#!/usr/bin/python
# M Orsi Oct 2009
# This script implements a free-energy based permeability predictor. The input file
# "free_energy.dat", normally obtained from z-constraint simulations, is processed to
# output the corresponding permeability prediction.
#
# Usage:
# python /home/orsi/brahmsProject/analysisCode/zConstraint/permPredictor_2sets.py G12.dat

import sys, string, math

try:
    file = sys.argv[1]
except IndexError:
    raise "Provide input file!"

inFile = open(file,'r')
lines=inFile.readlines()
inFile.close()

kB = 1.3806505e-26 # [kJ/K]
AVOGAD = 6.0221367e23 # [mol^(-1)]
kBT = kB * 303
kBT_mol = kBT*AVOGAD
#print kBT_mol

nSets = 2;
print "nSets = %d" % nSets

p = []

for i in range(nSets):
    maxDeltaG = minDeltaG = 0.; # initialise transfer free energy values
    for line in lines:
        words = string.split( line )
        deltaG = float(words[i+1]) # transfer free energy
        if deltaG > maxDeltaG:
            maxDeltaG = deltaG
        if deltaG < minDeltaG:
            minDeltaG = deltaG
    p.append(kBT_mol/(maxDeltaG-minDeltaG))

# p now contains "nSets" data points - now get average and standard error!

print p

sum = sum2 = sumSq = 0.
for i in range(nSets):
    #print i, p[i]
    sum = sum + p[i]
    sum2 = sum2 + p[i]**2
sum2 = sum2 / nSets
pAvg = sum / nSets
sumSq = pAvg**2
sd = (sum2 - sumSq)**0.5
se = sd/(nSets**0.5)

print "Permeability predictor = %.5f %.5f" % (pAvg,se)

#p = 1/(maxDeltaG-minDeltaG)
#p = kBT_mol/(math.e**(maxDeltaG-minDeltaG)) # NO! wrong units!!!
#p = 100*kBT_mol/(maxDeltaG-minDeltaG)
#p = math.e**(kBT_mol/(maxDeltaG-minDeltaG))
#p = 10**(kBT_mol/(maxDeltaG-minDeltaG))

#print "Permeability predictor: %.4g" % p
